clear

*Note: first change working directory replication package folder

run "scripts/programs/program_data.do"
run "scripts/programs//program_regs.do"
run "scripts/programs/program_regtable.do"


/*********************************************************************************/
*Set program parameters
local radius = 20
local angle = 15
local length = 1
/**********************************************************************************/

/*************************	DATA FOR SUMMARY STATISTICS ******************************/

data anyhu hu hudens ///
	anypop pop lnpop popdens inc_pcincome huxpcinc huxpcinc_p race_pct_white ///
	any_res count_res count_res_dens avgval_res val_res ///
	watershedimp hablinepct habpolypct ///
	roadpct majorroad tri ///
	tmissing tdiff lnt lnint campground wilderness wind_diff ///
	using "data/20-by-1km_15deg/data_centroid.dta", ///
	burnt(centroid) nodropout

merge m:1 firenum using ///
	"data/20-by-1km_15deg/sample_census.dta", ///
	keepusing(sample) nogen

/**********************************************************************************/
*Label variables

keep if sample == 1

gen _int = cond(lnint!=0, exp(lnint) - 1,.)

label var pop "Total population"
label var anyhu "(No. HU $> 0$)"
label var hudens "HU density (HU/sq. km)"
label var inc_pcincome "Per cap. income (thousands USD)"
label var huxpcinc "No. HU $\times$ per cap. inc. (millions USD)"
label var huxpcinc_p "Total income/sq. km"
label var race_pct_white "Pct. white"
label var watershedimp "Watershed importance"
label var roadpct "Pct. $<0.5 km$ fromm road"
label var majorroad "Contains major road"

label var toa "$ T $ (hours since ignition)"
label var tri "Topographic ruggedness index"
label var tdiff "$\Delta T$ (hours until arrival in next cell)"
label var tmissing "(No fuels)"
label var _int "Fire intensity (kW/hour)"
label var hablinepct "Pct. TES habitat (non-stream)"
label var habpolypct "TES habitat"
label var wilderness "Wilderness"
label var campground "Campgrounds"
label var wind_diff "Wind difference"
egen distquart = cut(distnum), at(0,6,11,16,21)
tab distquart, gen(dq)

/********************************************************************************/
/* PANEL 1 */
/********************************************************************************/

local varlist anyhu hu hudens inc_pcincome huxpcinc ///
	watershedimp habpolypct wilderness campground ///
	tri roadpct  ///
	toa tdiff tmissing _int majorroad wind_diff

replace inc_pcinc = . if inc_pcinc == 0

forvalues i = 1/4 {
	eststo dq`i': qui estpost summ `varlist' if dq`i' == 1 & sample == 1
	}

eststo all: estpost summ `varlist' if sample == 1

local output = "Results/Descriptives/table_contents.tex"

estout dq1 dq2 dq3 dq4 all using `output', replace style(tex) ///
	numbers("{\bf (" ")}") cells("mean(fmt(%12.3gc))") ///
	stats(N, label("\doubleT Number of obs.")) label ///
	mlabels("{\bf 0-5 km}" "{\bf 5-10 km}" "{\bf 10-15 km}" "{\bf 15-20 km}" "{\bf Whole Sample}", end("\midrule \T\B")) ///
	refcat(anyhu "\doubleT {\bf I. Full sample} &&&&& \\ \doubleT{\em Benefit vars.}"  ///
	toa "\T {\em Fire spread vars.}" ///
	tri "\T {\em Cost vars.}" ///
	watershedimp "\T {\em Other values at risk}" ///
	, nolabel) ///
	collabels(none)

/********************************************************************************/
/* PANEL 2 */
/********************************************************************************/

data anyhu hu hudens ///
	anypop pop lnpop popdens inc_pcincome huxpcinc huxpcinc_p race_pct_white ///
	any_res count_res count_res_dens avgval_res val_res ///
	watershedimp hablinepct habpolypct ///
	roadpct majorroad tri ///
	tmissing tdiff lnt lnint campground wilderness wind_diff ///
	using "data/20-by-1km_15deg/data_centroid.dta", ///
	burnt(centroid) nodropout

merge m:1 firenum using ///
	"data/20-by-1km_15deg/sample_assessors.dta", ///
	keepusing(sample) nogen

egen distquart = cut(distnum), at(0,6,11,16,21)
tab distquart, gen(dq)

label var any_res "\doubleT (No. res. props. $> 0$)"
label var count_res "No. residential props."
label var avgval_res "Avg. value res. props. (millions USD)"
label var val_res "Total value res. props. (millions USD)"

local varlist any_res count_res count_res_dens avgval_res val_res
replace avgval_res = . if count_res == 0

forvalues i = 1/4 {
	eststo dq`i': qui estpost summ `varlist' if dq`i' == 1 & sample == 1
	}

eststo all: estpost summ `varlist' if sample == 1

local output = "Results/Descriptives/panel_2.tex"

estout dq1 dq2 dq3 dq4 all using `output', replace style(tex) ///
	cells("mean(fmt(%12.3gc))") ///
	stats(N, label("\doubleT Number of obs.")) label ///
	mlabels(, none) ///
	refcat(any_res "\doubleT {\bf II. 2011-2015 Assessor's Data Sample}" , nolabel) ///
	collabels(none)
